# Open the Stata file
memory.limit(4000)
library(foreign)
library(Zelig)

## Set work drive to read in data file
#setwd("N:/Govt 2001/Hunt")

hunt.mat <- read.dta("Hunt_easyload.dta")

###########################################
#############CLEANING DATA#################
###########################################

## SUBSETTING BASED ON TIME##
## Drop obs when question wasn't asked (p 397 Hunt article)

hunt.mat <- hunt.mat[hunt.mat$year %in% c(1977, 1985, 1986, 1988, 1989, 1990, 1991, 1993, 1994, 1996, 1998, 2000, 2002, 2004, 2006, 2008, 2010), ]

###############################
##CREATING DEPENDENT VARIABLE##
###############################

## Rename racdif1-4 for clarity - Discrimination is racdif1, ability is racdif2, education is racdif3, motivation is racdif4

library(gregmisc)
hunt.mat <- rename.vars(hunt.mat, c("racdif1", "racdif2","racdif3", "racdif4"), c("racdif_disc", "racdif_abil", "racdif_educ", "racdif_motiv"))

## Generate categories of explanations for BW gap
hunt.mat$ability <- as.numeric(hunt.mat$racdif_abil=="yes" & hunt.mat$racdif_educ == "no" & hunt.mat$racdif_disc == "no")
hunt.mat$motivation <- as.numeric(hunt.mat$racdif_motiv=="yes" & hunt.mat$racdif_educ == "no" & hunt.mat$racdif_disc == "no" & hunt.mat$racdif_abil=="no")
hunt.mat$abil_struct <- as.numeric(hunt.mat$racdif_abil=="yes" & (hunt.mat$racdif_educ == "yes" | hunt.mat$racdif_disc == "yes"))
hunt.mat$motiv_struct <- as.numeric(hunt.mat$racdif_abil=="no" & hunt.mat$racdif_motiv=="yes" & (hunt.mat$racdif_educ == "yes" | hunt.mat$racdif_disc == "yes"))
hunt.mat$education <- as.numeric(hunt.mat$racdif_abil=="no" & hunt.mat$racdif_motiv=="no" & hunt.mat$racdif_educ == "yes" & hunt.mat$racdif_disc == "no")
hunt.mat$discrimination <- as.numeric(hunt.mat$racdif_abil=="no" & hunt.mat$racdif_motiv=="no" & hunt.mat$racdif_disc == "yes")
hunt.mat$none <- as.numeric(hunt.mat$racdif_abil=="no" & hunt.mat$racdif_motiv=="no" & hunt.mat$racdif_educ == "no" & hunt.mat$racdif_disc == "no")


## Turn the racdif variables into 0/1 dummies (useful for multiple imputation later)
hunt.mat$racdif_disc <- as.numeric(hunt.mat$racdif_disc=="yes")
hunt.mat$racdif_abil <- as.numeric(hunt.mat$racdif_abil=="yes")
hunt.mat$racdif_educ <- as.numeric(hunt.mat$racdif_educ=="yes")
hunt.mat$racdif_motiv <- as.numeric(hunt.mat$racdif_motiv == "yes")

##################################
##CREATING INDEPENDENT VARIABLES##
##################################

###########RACE VARIABLE##########
hunt.mat$hunt_race <- hunt.mat$race		#assigning race to hunt_race because race already captures those who are white or black

## Because hunt_race is a factor, I need to add "hispanic" to the levels otherwise it won't accept it
hunt.mat$hunt_race <- factor(hunt.mat$hunt_race, levels = c(levels(hunt.mat$hunt_race), "hispanic"))

hunt.mat$hunt_race[hunt.mat$ethnic %in% c(17, 22, 25, 38)] <- "hispanic"
hunt.mat$hunt_race[hunt.mat$ethnum=="cannot choose 1" & (hunt.mat$eth1 %in% c(17, 22, 25, 38) & hunt.mat$eth2 %in% c(17, 22, 25, 38))] <- "hispanic"
hunt.mat$hunt_race[hunt.mat$ethnum=="cannot choose 1" & (hunt.mat$eth1 %in% c(17, 22, 25, 38)& hunt.mat$eth3 %in% c(17, 22, 25, 38))] <- "hispanic"
hunt.mat$hunt_race[hunt.mat$ethnum=="cannot choose 1" & (hunt.mat$eth3 %in% c(17, 22, 25, 38) & hunt.mat$eth2 %in% c(17, 22, 25, 38))] <- "hispanic"

table(hunt.mat$hunt_race, exclude=NULL)

## Remove obs for race=other
hunt.mat <- hunt.mat[!hunt.mat$hunt_race=="other",]

## Eliminate category "iap": first assign its members to missing, then drop the category itself
is.na(hunt.mat$hunt_race)<- hunt.mat$hunt_race=="iap"
hunt.mat$hunt_race <- factor(hunt.mat$hunt_race, levels = c("white","black","hispanic"))

## Generate dummy variables
hunt.mat$black <- as.numeric(hunt.mat$hunt_race=="black")
hunt.mat$hispanic <- as.numeric(hunt.mat$hunt_race=="hispanic")
hunt.mat$white <- as.numeric(hunt.mat$hunt_race=="white")

###########CLASS VARIABLE##########
## Check the variable
table(hunt.mat$class, exclude=NULL)

## Let's reclassify the 1 "no class" as NA and get rid of the empty levels
hunt.mat$hunt_class <- hunt.mat$class
is.na(hunt.mat$hunt_class)<- hunt.mat$hunt_class=="no class"
hunt.mat$hunt_class <- factor(hunt.mat$hunt_class, levels = c("lower class","working class","middle class", "upper class"))

## Check
table(hunt.mat$hunt_class, exclude=NULL) #OK

## Generate dummy variables
hunt.mat$lowclass <- as.numeric(hunt.mat$class=="lower class")
hunt.mat$workclass<- as.numeric(hunt.mat$class=="working class")

###########AGE VARIABLE##########
## Age is already generated (age)
 
###########REGION VARIABLE##########
## Generate region where south is south atlantic, e. south central or w. south central or zero otherwise (no missing in region initial variable)
hunt.mat$hunt_region <- as.numeric(hunt.mat$region == "south atlantic" | hunt.mat$region=="e. sou. central" | hunt.mat$region == "w. sou. central")
table(hunt.mat$hunt_region)

###########POLITICAL VIEWS VARIABLE##########
table(hunt.mat$polviews, exclude=NULL)

## Recoding as numbers because R only recognizes as text
hunt.mat$hunt_polviews <- as.numeric(hunt.mat$polviews)
table(hunt.mat$hunt_polviews, exclude=NULL)

## Drops iap, dk, na and numbers substantive categories 2 to 8. Also outputs 4076 missing so no pb there.
#Now change to a scale of 1 to 7 
hunt.mat$hunt_polviews <- hunt.mat$hunt_polviews - 1

table(hunt.mat$hunt_polviews, exclude=NULL) #OK


###########SEX DUMMY VARIABLE##########
hunt.mat$female <- as.numeric(hunt.mat$sex =="female")


###########FUNDAMENTALIST VARIABLE##########
## Generate fundamentalist variable where it's 1 if person is a fundamentalist, 0  if otherwise (moderate or liberal)
hunt.mat$hunt_fund <- as.numeric(hunt.mat$fund == "fundamentalist")
table(hunt.mat$hunt_fund, exclude=NULL)


###########YEAR VARIABLE##########
## Year variable is already generated (year) - no missing values


###########EDUCATION VARIABLE##########
## Education variable is already generated (educ)


##########################################################
#Reconstructing outcome variables into �person-centered,�# 
#�mixed� or �structuralist.�                             #
##########################################################

## Person centered = person_cent
## Mixed modes = mixed
## structuralist = struct

hunt.mat$person_cent <- as.numeric(hunt.mat$ability == 1 | hunt.mat$motivation == 1)
hunt.mat$mixed <- as.numeric(hunt.mat$abil_struct == 1 | hunt.mat$motiv_struct == 1)
hunt.mat$struct <- as.numeric(hunt.mat$education == 1 | hunt.mat$discrimination == 1)

##############SEPARATING HUNT AND OBAMA DATASETS########################
## Careful, the order is important
obama.mat <- hunt.mat
hunt.mat <- hunt.mat[hunt.mat$year %in% c(1977, 1985, 1986, 1988, 1989, 1990, 1991, 1993, 1994, 1996, 1998, 2000, 2002, 2004), ]

## Create a dummy for post-2008 in Obama dataset
obama.mat$obama <- as.numeric(obama.mat$year>=2008)


###########################################
########## Replicating Table 4 ############
###########################################

###########Listwise deletion##########
#Keep only the variables we'll use in the analysis and remove all missing data
hunt.listwise<- subset(hunt.mat[hunt.mat$year!=1977,], select = c(hunt_race, ability, motivation, 
	abil_struct, motiv_struct, education, discrimination, none, person_cent, mixed, struct, black, 
	hispanic,educ, lowclass, workclass, age, female, hunt_region, hunt_fund, 
	hunt_polviews, year, wtss))
hunt.listwise <- na.omit(hunt.listwise)
dim(hunt.listwise)


########Logit with glm#########

fit.ability <- glm(ability~black+hispanic+educ+lowclass+workclass+age+female+hunt_region+hunt_fund+hunt_polviews+year, family = binomial(link="logit"), data=hunt.listwise, weights=wtss)
fit.motivation <- glm(motivation~black+hispanic+educ+lowclass+workclass+age+female+hunt_region+hunt_fund+hunt_polviews+year, family = binomial(link="logit"), data=hunt.listwise, weights=wtss)
fit.abil_struct <- glm(abil_struct~black+hispanic+educ+lowclass+workclass+age+female+hunt_region+hunt_fund+hunt_polviews+year, family = binomial(link="logit"), data=hunt.listwise, weights=wtss)
fit.motiv_struct <- glm(motiv_struct~black+hispanic+educ+lowclass+workclass+age+female+hunt_region+hunt_fund+hunt_polviews+year, family = binomial(link="logit"), data=hunt.listwise, weights=wtss)
fit.education <- glm(education~black+hispanic+educ+lowclass+workclass+age+female+hunt_region+hunt_fund+hunt_polviews+year, family = binomial(link="logit"), data=hunt.listwise, weights=wtss)
fit.discrimination <- glm(discrimination~black+hispanic+educ+lowclass+workclass+age+female+hunt_region+hunt_fund+hunt_polviews+year, family = binomial(link="logit"), data=hunt.listwise, weights=wtss)
fit.none <- glm(none~black+hispanic+educ+lowclass+workclass+age+female+hunt_region+hunt_fund+hunt_polviews+year, family = binomial(link="logit"), data=hunt.listwise, weights=wtss)

########Creating standard errors##########
## Install.packages("sandwich")
library(sandwich) 	

## Calculating standard errors for each model generated above
bread.ability <- vcov(fit.ability)
est.fun.ability <- estfun(fit.ability)		#create matrices first derivative
meat.ability <- t(est.fun.ability)%*% est.fun.ability
sandwich.ability <- bread.ability %*% meat.ability %*% bread.ability
## Outputs robust vcov matrix

bread.motivation <- vcov(fit.motivation)
est.fun.motivation <- estfun(fit.motivation)		
meat.motivation <- t(est.fun.motivation)%*% est.fun.motivation
sandwich.motivation <- bread.motivation %*% meat.motivation %*% bread.motivation

bread.abil_struct <- vcov(fit.abil_struct)
est.fun.abil_struct<- estfun(fit.abil_struct)		
meat.abil_struct<- t(est.fun.abil_struct)%*% est.fun.abil_struct
sandwich.abil_struct<- bread.motivation %*% meat.abil_struct%*% bread.abil_struct

bread.motiv_struct <- vcov(fit.motiv_struct)
est.fun.motiv_struct<- estfun(fit.motiv_struct)		
meat.motiv_struct<- t(est.fun.motiv_struct)%*% est.fun.motiv_struct
sandwich.motiv_struct<- bread.motiv_struct%*% meat.motiv_struct%*% bread.motiv_struct

bread.education <- vcov(fit.education )
est.fun.education <- estfun(fit.education )		
meat.education <- t(est.fun.education )%*% est.fun.education 
sandwich.education <- bread.education %*% meat.education %*% bread.education 

bread.discrimination <- vcov(fit.discrimination )
est.fun.discrimination <- estfun(fit.discrimination )		
meat.discrimination <- t(est.fun.discrimination )%*% est.fun.discrimination 
sandwich.discrimination <- bread.discrimination %*% meat.discrimination %*% bread.discrimination 

bread.none <- vcov(fit.none )
est.fun.none <- estfun(fit.none )		
meat.none <- t(est.fun.none )%*% est.fun.none 
sandwich.none <- bread.none %*% meat.none %*% bread.none 

#############Creating p-values#############
## Install.packages("lmtest")
library(lmtest)
coef.ability <- coeftest(fit.ability, sandwich.ability)
coef.motivation <- coeftest(fit.motivation, sandwich.motivation)
coef.abil_struct <- coeftest(fit.abil_struct, sandwich.abil_struct)
coef.motiv_struct <- coeftest(fit.motiv_struct, sandwich.motiv_struct)
coef.education <- coeftest(fit.education, sandwich.education)
coef.discrimination <- coeftest(fit.discrimination, sandwich.discrimination)
coef.none <- coeftest(fit.none, sandwich.none)

############Outputting Table 4##############

library(xtable)

## Storing results from p-values
tab.ability <- coef.ability[,1:2]
tab.motivation <- coef.motivation[,1:2]
tab.abil_struct <- coef.abil_struct[,1:2]
tab.motiv_struct <- coef.motiv_struct[,1:2]
tab.education <- coef.education[,1:2]
tab.discrimination <- coef.discrimination[,1:2]
tab.none <- coef.none[,1:2]

table4 <- cbind(tab.ability, tab.motivation, tab.abil_struct, tab.motiv_struct, tab.education, tab.discrimination, tab.none)
rows <- nrow(table4)
table4 <- table4[2:rows,]			#dropping intercept

View(table4)

## Write table to excel and format table
table4 <- round(table4, digits=3)
rownames(table4) <-c("Black", "Hispanic", "Education", "Lower Class", "Working Class", "Age", 
	"Female", "South", "Fundamentalist", "Conservative", "Year")

write.csv(table4, file="Table_4.csv",row.names=T)

#######################################################
######EXPLORING POTENTIAL IMPROVEMENTS TO HUNT 07######
#######################################################

########################################################
#Running analysis on aggregate modes �person-centered,�# 
#�mixed� or �structuralist.�                           #
########################################################

## Running logit for dummies
fit.person_cent <- glm(person_cent~black+hispanic+educ+lowclass+workclass+age+female+hunt_region+hunt_fund+hunt_polviews+year, family = binomial(link="logit"), data=hunt.listwise, weights=wtss)
fit.mixed <- glm(mixed~black+hispanic+educ+lowclass+workclass+age+female+hunt_region+hunt_fund+hunt_polviews+year, family = binomial(link="logit"), data=hunt.listwise, weights=wtss)
fit.struct <- glm(struct~black+hispanic+educ+lowclass+workclass+age+female+hunt_region+hunt_fund+hunt_polviews+year, family = binomial(link="logit"), data=hunt.listwise, weights=wtss)

summary(fit.person_cent)
summary(fit.mixed)
summary(fit.struct)


#######################################################
#############Quantities of Interest####################
#######################################################

## We present them over time, so define:
year.hunt <- c(1985,1986, 1988, 1989,1990, 1991, 1993, 1994,1996, 1998, 2000, 2002, 2004)

## Set predictive covariates to median - model used odes not matter since they have the same covariates
	# Allow to vary by year
Xwhite.hunt <-  setx(fit.person_cent, fn=list(numeric=median), black  =  0, hispanic=0, year=year.hunt)
Xblack.hunt <-  setx(fit.person_cent, fn=list(numeric=median),  black  =  1, hispanic=0, year=year.hunt)
Xhisp.hunt  <-  setx(fit.person_cent, fn=list(numeric=median),  black  =  0, hispanic=1, year=year.hunt)

## Since we cannot use sim when models are not run with Zelig, compute first difference manually
	#Turn our X's into matrices
Xwhite.hunt <- as.matrix(Xwhite.hunt)
Xblack.hunt <- as.matrix(Xblack.hunt)
Xhisp.hunt <- as.matrix(Xhisp.hunt)


#################
#Person-centered#
#################

## Simulate betas from logit results
set.seed(12345)
beta.hunt.pcent <- mvrnorm(10000, mu = fit.person_cent$coef, Sigma =vcov(fit.person_cent))

## Calculate first differences for blacks v. whites
fd.hunt.pcentblack <- (1 / (1 + exp(-Xblack.hunt%*%t(beta.hunt.pcent)))) - (1 / (1 + exp(-Xwhite.hunt%*%t(beta.hunt.pcent))))
meanfd.hunt.pcentblack <- apply(fd.hunt.pcentblack, 1, mean)
cilowfd.hunt.pcentblack<- apply(fd.hunt.pcentblack, 1, quantile, .025)
cihighfd.hunt.pcentblack <- apply(fd.hunt.pcentblack, 1, quantile, .975)

## Calculate first differences for hisp v. whites
fd.hunt.pcenthisp<- (1 / (1 + exp(-Xhisp.hunt%*%t(beta.hunt.pcent)))) - (1 / (1 + exp(-Xwhite.hunt%*%t(beta.hunt.pcent))))
meanfd.hunt.pcenthisp <- apply(fd.hunt.pcenthisp, 1, mean)
cilowfd.hunt.pcenthisp <- apply(fd.hunt.pcenthisp, 1, quantile, .025)
cihighfd.hunt.pcenthisp <- apply(fd.hunt.pcenthisp, 1, quantile, .975)

## Plot with vertical CI for each year - Person-centered
year <- year.hunt
plot(x=year ,  xlab="Year", y=meanfd.hunt.pcentblack, , col  =  "blue", main="", sub="", ylab="First differences", ylim  =  c(-0.25,0.01))
title(main = list("Figure 1: Effect of race on probability of person-centered explanation (1985-2004)", cex=1.25, font=2), sub=list("Original specification", cex=1, font=2))
segments(x0= year, x1= year, y0  =  cilowfd.hunt.pcentblack, y1  =  cihighfd.hunt.pcentblack, col  =  "blue")
points(x=year ,  y=meanfd.hunt.pcenthisp,  col  =  "red")
segments(x0= year, x1= year, y0  =  cilowfd.hunt.pcenthisp, y1  =  cihighfd.hunt.pcenthisp, col  =  "red")

legend("topright", c("Blacks v. Whites","Hispanics v. Whites"),
  col=c("blue", "red"), lty=c(1,1))

#######
#Mixed#
#######

## Simulate betas from logit results
set.seed(12345)
beta.hunt.mixed <- mvrnorm(10000, mu = fit.mixed$coef, Sigma =vcov(fit.mixed))

## Calculate first differences for blacks v. whites
fd.hunt.mixedblack <- (1 / (1 + exp(-Xblack.hunt%*%t(beta.hunt.mixed)))) - (1 / (1 + exp(-Xwhite.hunt%*%t(beta.hunt.mixed))))
meanfd.hunt.mixedblack <- apply(fd.hunt.mixedblack, 1, mean)
cilowfd.hunt.mixedblack<- apply(fd.hunt.mixedblack, 1, quantile, .025)
cihighfd.hunt.mixedblack <- apply(fd.hunt.mixedblack, 1, quantile, .975)

## Calculate first differences for hisp v. whites
fd.hunt.mixedhisp<- (1 / (1 + exp(-Xhisp.hunt%*%t(beta.hunt.mixed)))) - (1 / (1 + exp(-Xwhite.hunt%*%t(beta.hunt.mixed))))
meanfd.hunt.mixedhisp <- apply(fd.hunt.mixedhisp, 1, mean)
cilowfd.hunt.mixedhisp <- apply(fd.hunt.mixedhisp, 1, quantile, .025)
cihighfd.hunt.mixedhisp <- apply(fd.hunt.mixedhisp, 1, quantile, .975)

## Plot with vertical CI for each year - Mixed
plot(x=year ,  xlab="Year", y=meanfd.hunt.mixedblack, , col  =  "blue", main="", sub="", ylab="First differences", ylim  =  c(-0.05,0.25))
title(main = list("Figure 1: Effect of race on probability of mixed explanation (1985-2004)", cex=1.25, font=2), sub=list("Original specification", cex=1, font=2))
segments(x0= year, x1= year, y0  =  cilowfd.hunt.mixedblack, y1  =  cihighfd.hunt.mixedblack, col  =  "blue")
points(x=year ,  y=meanfd.hunt.mixedhisp,  col  =  "red")
segments(x0= year, x1= year, y0  =  cilowfd.hunt.mixedhisp, y1  =  cihighfd.hunt.mixedhisp, col  =  "red")

legend("topright", c("Blacks v. Whites","Hispanics v. Whites"),
  col=c("blue", "red"), lty=c(1,1))

############
#Structural#
############

## Simulate betas from logit results
set.seed(12345)
beta.hunt.struct <- mvrnorm(10000, mu = fit.struct$coef, Sigma =vcov(fit.struct))

## Calculate first differences for blacks v. whites
fd.hunt.structblack <- (1 / (1 + exp(-Xblack.hunt%*%t(beta.hunt.struct)))) - (1 / (1 + exp(-Xwhite.hunt%*%t(beta.hunt.struct))))
meanfd.hunt.structblack <- apply(fd.hunt.structblack, 1, mean)
cilowfd.hunt.structblack<- apply(fd.hunt.structblack, 1, quantile, .025)
cihighfd.hunt.structblack <- apply(fd.hunt.structblack, 1, quantile, .975)

## Calculate first differences for hisp v. whites
fd.hunt.structhisp<- (1 / (1 + exp(-Xhisp.hunt%*%t(beta.hunt.struct)))) - (1 / (1 + exp(-Xwhite.hunt%*%t(beta.hunt.struct))))
meanfd.hunt.structhisp <- apply(fd.hunt.structhisp, 1, mean)
cilowfd.hunt.structhisp <- apply(fd.hunt.structhisp, 1, quantile, .025)
cihighfd.hunt.structhisp <- apply(fd.hunt.structhisp, 1, quantile, .975)


## Plot with vertical CI for each year - Structural
plot(x=year ,  xlab="Year", y=meanfd.hunt.structblack, , col  =  "blue", main="", sub="", ylab="First differences", ylim  =  c(-0.07,0.4))
title(main = list("Figure 1: Effect of race on probability of structural explanation (1985-2004)", cex=1.25, font=2), sub=list("Original specification", cex=1, font=2))
segments(x0= year, x1= year, y0  =  cilowfd.hunt.structblack, y1  =  cihighfd.hunt.structblack, col  =  "blue")
points(x=year ,  y=meanfd.hunt.structhisp,  col  =  "red")
segments(x0= year, x1= year, y0  =  cilowfd.hunt.structhisp, y1  =  cihighfd.hunt.structhisp, col  =  "red")

legend("topright", c("Blacks v. Whites","Hispanics v. Whites"),
  col=c("blue", "red"), lty=c(1,1))


#########################################
########Multiple imputation##############
#########################################
set.seed(12345)
require(Amelia)
library(Zelig)

## Creating dataset for amelia that only has variables we include in the regression
amelia.data <- with(hunt.mat, as.data.frame(cbind(racdif_abil, racdif_motiv, racdif_educ, 
	racdif_disc, black, hispanic,educ, lowclass, workclass, age, female, hunt_region, 
	hunt_fund, hunt_polviews, year)))
amelia.data <- amelia.data[amelia.data$year!=1977,]

## noms lets us set which variables are categorical
noms <- c("racdif_disc", "racdif_abil", "racdif_educ", "racdif_motiv", "hunt_region", "hunt_fund", "hunt_polviews")

## We use the amelia() package to create the five bootstrapped em imputation models.  
	# amelia.out will be a list of 5 imputed datasets, each of which can be accessed using 
	# a.out$imputations[[i]]
amelia.out <- amelia(x = amelia.data, m = 5, noms = noms) # noms tells r which variables are nominal
summary(amelia.out$imputations)
View(amelia.out$imputations$imp1)
dim(amelia.out$imputations$imp1)

## Save imputed datasets (default .csv)
write.amelia(obj=amelia.out, file.stem = "amelia.out")

## Following section covers reading in files so we don't have to run imputations each time. Do not run
	## if generating new imputations
amelia.out <- c()
amelia.out$imputations <- c()
amelia.out$imputations$imp1 <- read.csv("amelia.out1.csv")
amelia.out$imputations$imp2 <- read.csv("amelia.out2.csv")
amelia.out$imputations$imp3 <- read.csv("amelia.out3.csv")
amelia.out$imputations$imp4 <- read.csv("amelia.out4.csv")
amelia.out$imputations$imp5 <- read.csv("amelia.out5.csv")

## Getting rid of first column, X, which is just the row number
amelia.out$imputations$imp1 <- subset(amelia.out$imputations$imp1, select = c(racdif_abil, racdif_motiv, racdif_educ, 
	racdif_disc, black, hispanic, educ, lowclass, workclass, age, female, hunt_region, hunt_fund, hunt_polviews, year))
amelia.out$imputations$imp2 <- subset(amelia.out$imputations$imp2, select = c(racdif_abil, racdif_motiv, racdif_educ, 
	racdif_disc, black, hispanic, educ, lowclass, workclass, age, female, hunt_region, hunt_fund, hunt_polviews, year))
amelia.out$imputations$imp3 <- subset(amelia.out$imputations$imp3, select = c(racdif_abil, racdif_motiv, racdif_educ, 
	racdif_disc, black, hispanic, educ, lowclass, workclass, age, female, hunt_region, hunt_fund, hunt_polviews, year))
amelia.out$imputations$imp4 <- subset(amelia.out$imputations$imp4, select = c(racdif_abil, racdif_motiv, racdif_educ, 
	racdif_disc, black, hispanic, educ, lowclass, workclass, age, female, hunt_region, hunt_fund, hunt_polviews, year))
amelia.out$imputations$imp5 <- subset(amelia.out$imputations$imp5, select = c(racdif_abil, racdif_motiv, racdif_educ, 
	racdif_disc, black, hispanic, educ, lowclass, workclass, age, female, hunt_region, hunt_fund, hunt_polviews, year))
## End of section for reading in previously imputed datasets 

## Defining ability
amelia.out$imputations$imp1$ability <- with(amelia.out$imputations$imp1, as.numeric(racdif_abil==1 & racdif_educ == 0 & racdif_disc == 0))
amelia.out$imputations$imp2$ability <- with(amelia.out$imputations$imp2, as.numeric(racdif_abil==1 & racdif_educ == 0 & racdif_disc == 0))
amelia.out$imputations$imp3$ability <- with(amelia.out$imputations$imp3, as.numeric(racdif_abil==1 & racdif_educ == 0 & racdif_disc == 0))
amelia.out$imputations$imp4$ability <- with(amelia.out$imputations$imp4, as.numeric(racdif_abil==1 & racdif_educ == 0 & racdif_disc == 0))
amelia.out$imputations$imp5$ability <- with(amelia.out$imputations$imp5, as.numeric(racdif_abil==1 & racdif_educ == 0 & racdif_disc == 0))

## Defining motivation
amelia.out$imputations$imp1$motivation <- with(amelia.out$imputations$imp1, as.numeric(racdif_motiv== 1 & racdif_educ == 0 & racdif_disc == 0 & racdif_abil==0))
amelia.out$imputations$imp2$motivation <- with(amelia.out$imputations$imp2, as.numeric(racdif_motiv== 1 & racdif_educ == 0 & racdif_disc == 0 & racdif_abil==0))
amelia.out$imputations$imp3$motivation <- with(amelia.out$imputations$imp3, as.numeric(racdif_motiv== 1 & racdif_educ == 0 & racdif_disc == 0 & racdif_abil==0))
amelia.out$imputations$imp4$motivation <- with(amelia.out$imputations$imp4, as.numeric(racdif_motiv== 1 & racdif_educ == 0 & racdif_disc == 0 & racdif_abil==0))
amelia.out$imputations$imp5$motivation <- with(amelia.out$imputations$imp5, as.numeric(racdif_motiv== 1 & racdif_educ == 0 & racdif_disc == 0 & racdif_abil==0))

## Defining abil_struct
amelia.out$imputations$imp1$abil_struct <- with(amelia.out$imputations$imp1, as.numeric(racdif_abil== 1 & (racdif_educ == 1 | racdif_disc == 1)))
amelia.out$imputations$imp2$abil_struct <- with(amelia.out$imputations$imp2, as.numeric(racdif_abil== 1 & (racdif_educ == 1 | racdif_disc == 1)))
amelia.out$imputations$imp3$abil_struct <- with(amelia.out$imputations$imp3, as.numeric(racdif_abil== 1 & (racdif_educ == 1 | racdif_disc == 1)))
amelia.out$imputations$imp4$abil_struct <- with(amelia.out$imputations$imp4, as.numeric(racdif_abil== 1 & (racdif_educ == 1 | racdif_disc == 1)))	
amelia.out$imputations$imp5$abil_struct <- with(amelia.out$imputations$imp5, as.numeric(racdif_abil== 1 & (racdif_educ == 1 | racdif_disc == 1)))

## Defining motiv_struct
amelia.out$imputations$imp1$motiv_struct <- with(amelia.out$imputations$imp1, as.numeric(racdif_abil== 0 & racdif_motiv== 1 & (racdif_educ == 1 | racdif_disc == 1)))
amelia.out$imputations$imp2$motiv_struct <- with(amelia.out$imputations$imp2, as.numeric(racdif_abil== 0 & racdif_motiv== 1 & (racdif_educ == 1 | racdif_disc == 1)))
amelia.out$imputations$imp3$motiv_struct <- with(amelia.out$imputations$imp3, as.numeric(racdif_abil== 0 & racdif_motiv== 1 & (racdif_educ == 1 | racdif_disc == 1)))
amelia.out$imputations$imp4$motiv_struct <- with(amelia.out$imputations$imp4, as.numeric(racdif_abil== 0 & racdif_motiv== 1 & (racdif_educ == 1 | racdif_disc == 1)))
amelia.out$imputations$imp5$motiv_struct <- with(amelia.out$imputations$imp5, as.numeric(racdif_abil== 0 & racdif_motiv== 1 & (racdif_educ == 1 | racdif_disc == 1)))

## Defining education
amelia.out$imputations$imp1$education <- with(amelia.out$imputations$imp1, as.numeric(racdif_abil== 0 & racdif_motiv==0 & racdif_educ == 1 & racdif_disc == 0))
amelia.out$imputations$imp2$education <- with(amelia.out$imputations$imp2, as.numeric(racdif_abil== 0 & racdif_motiv==0 & racdif_educ == 1 & racdif_disc == 0))
amelia.out$imputations$imp3$education <- with(amelia.out$imputations$imp3, as.numeric(racdif_abil== 0 & racdif_motiv==0 & racdif_educ == 1 & racdif_disc == 0))
amelia.out$imputations$imp4$education <- with(amelia.out$imputations$imp4, as.numeric(racdif_abil== 0 & racdif_motiv==0 & racdif_educ == 1 & racdif_disc == 0))
amelia.out$imputations$imp5$education <- with(amelia.out$imputations$imp5, as.numeric(racdif_abil== 0 & racdif_motiv==0 & racdif_educ == 1 & racdif_disc == 0))

## Defining discrimination
amelia.out$imputations$imp1$discrimination <- with(amelia.out$imputations$imp1, as.numeric(racdif_abil== 0 & racdif_motiv== 0 & racdif_disc == 1))
amelia.out$imputations$imp2$discrimination <- with(amelia.out$imputations$imp2, as.numeric(racdif_abil== 0 & racdif_motiv== 0 & racdif_disc == 1))
amelia.out$imputations$imp3$discrimination <- with(amelia.out$imputations$imp3, as.numeric(racdif_abil== 0 & racdif_motiv== 0 & racdif_disc == 1))
amelia.out$imputations$imp4$discrimination <- with(amelia.out$imputations$imp4, as.numeric(racdif_abil== 0 & racdif_motiv== 0 & racdif_disc == 1))
amelia.out$imputations$imp5$discrimination <- with(amelia.out$imputations$imp5, as.numeric(racdif_abil== 0 & racdif_motiv== 0 & racdif_disc == 1))

## Defining none
amelia.out$imputations$imp1$none <- with(amelia.out$imputations$imp1, as.numeric(racdif_abil == 0 & racdif_motiv == 0 & racdif_educ == 0 & racdif_disc == 0))
amelia.out$imputations$imp2$none <- with(amelia.out$imputations$imp2, as.numeric(racdif_abil == 0 & racdif_motiv == 0 & racdif_educ == 0 & racdif_disc == 0))
amelia.out$imputations$imp3$none <- with(amelia.out$imputations$imp3, as.numeric(racdif_abil == 0 & racdif_motiv == 0 & racdif_educ == 0 & racdif_disc == 0))
amelia.out$imputations$imp4$none <- with(amelia.out$imputations$imp4, as.numeric(racdif_abil == 0 & racdif_motiv == 0 & racdif_educ == 0 & racdif_disc == 0))
amelia.out$imputations$imp5$none <- with(amelia.out$imputations$imp5, as.numeric(racdif_abil == 0 & racdif_motiv == 0 & racdif_educ == 0 & racdif_disc == 0))


## Run models again, but using Zelig since it can handle multiply imputed data sets
imp.ability <- zelig(ability~black+hispanic+educ+lowclass+workclass+age+female+hunt_region+hunt_fund+hunt_polviews+year, family = binomial(link="logit"), data=amelia.out$imputations, model = "logit")
imp.motivation <- zelig(motivation~black+hispanic+educ+lowclass+workclass+age+female+hunt_region+hunt_fund+hunt_polviews+year, family = binomial(link="logit"), data=amelia.out$imputations, model = "logit")
imp.abil_struct <- zelig(abil_struct~black+hispanic+educ+lowclass+workclass+age+female+hunt_region+hunt_fund+hunt_polviews+year, family = binomial(link="logit"), data=amelia.out$imputations, model = "logit")
imp.motiv_struct <- zelig(motiv_struct~black+hispanic+educ+lowclass+workclass+age+female+hunt_region+hunt_fund+hunt_polviews+year, family = binomial(link="logit"), data=amelia.out$imputations, model = "logit")
imp.education <- zelig(education~black+hispanic+educ+lowclass+workclass+age+female+hunt_region+hunt_fund+hunt_polviews+year, family = binomial(link="logit"), data=amelia.out$imputations, model = "logit")
imp.discrimination <- zelig(discrimination~black+hispanic+educ+lowclass+workclass+age+female+hunt_region+hunt_fund+hunt_polviews+year, family = binomial(link="logit"), data=amelia.out$imputations, model = "logit")
imp.none <- zelig(none~black+hispanic+educ+lowclass+workclass+age+female+hunt_region+hunt_fund+hunt_polviews+year, family = binomial(link="logit"), data=amelia.out$imputations, model = "logit")


## The missingness map gives an overall sense of the shape and extent of the missingness.
missmap(amelia.out)

## The disperse function starts the algorithm at some unlikely values to check that 
	## AMELIA II hasn't found a local rather than a global maximum for the likelihood of 
	## the complete data
disperse(a.out, dims = 1, m = 5) 

## Generating table to compare to hunt's results for each outcome model
coef.imp.person_cent <- summary(imp.person_cent, subset = 1:5)$coefficients
coef.listwise.person_cent <- summary(fit.person_cent)$coefficients
table.person_cent <- write.csv(cbind(coef.imp.person_cent,"", coef.listwise.person_cent), 
	"Imputed versus listwise - Person-centered.csv",row.names=T)

coef.imp.mixed <- summary(imp.mixed, subset = 1:5)$coefficients
coef.listwise.mixed <- summary(fit.mixed)$coefficients
table.mixed<- write.csv(cbind(coef.imp.mixed,"", coef.listwise.mixed), 
	"Imputed versus listwise - Mixed.csv",row.names=T)

coef.imp.struct <- summary(imp.struct, subset = 1:5)$coefficients
coef.listwise.struct<- summary(fit.struct)$coefficients
table.struct<- write.csv(cbind(coef.imp.abil_struct,"", coef.listwise.abil_struct), 
	"Imputed versus listwise - Structural.csv",row.names=T)

################################################################
####################Obama Effect################################
################################################################

#######################################################
Reconstructing outcome variables into �person-centered,� 
�mixed� or �structuralist.� 
#######################################################
## Person centered = person_cent
## Mixed modes = mixed
## structuralist = struct

obama.mat$person_cent <- as.numeric(obama.mat$ability == 1 | obama.mat$motivation == 1)
obama.mat$mixed <- as.numeric(obama.mat$abil_struct == 1 | obama.mat$motiv_struct == 1)
obama.mat$struct <- as.numeric(obama.mat$education == 1 | obama.mat$discrimination == 1)


################################################
########## ANALYSIS for obama effect############
################################################

######Listwise deletion#########
obama.listwise <- subset(obama.mat[obama.mat$year!=1977,], select = c(person_cent, mixed, struct, obama, hunt_race, ability, motivation, 
	abil_struct, motiv_struct, education, discrimination, none, black, 
	hispanic,educ, lowclass, workclass, age, female, hunt_region, hunt_fund, 
	hunt_polviews, year, wtss))
obama.listwise <- na.omit(obama.listwise)
dim(obama.listwise)

########Running logit for 3 modes of explanation #########

## With the same specification as Hunt
## GLM binomial (probably underestimates Std. err but fits via log-likelihood so we can test fit)
obama.person_cent.same <- glm(person_cent~black+hispanic+educ+lowclass+workclass+age+female+hunt_region+hunt_fund+hunt_polviews+year, family = binomial(link="logit"), data=obama.listwise, weights=wtss)
obama.mixed.same <- glm(mixed~black+hispanic+educ+lowclass+workclass+age+female+hunt_region+hunt_fund+hunt_polviews+year, family = binomial(link="logit"), data=obama.listwise, weights=wtss)
obama.struct.same <- glm(struct~black+hispanic+educ+lowclass+workclass+age+female+hunt_region+hunt_fund+hunt_polviews+year, family = binomial(link="logit"), data=obama.listwise, weights=wtss)


## Run Obama specification
obama.person_cent <- glm(person_cent~black+hispanic+obama+(black*obama)+(hispanic*obama)+educ+lowclass+workclass+age+female+hunt_region+hunt_fund+hunt_polviews+year, family = binomial(link="logit"), data=obama.listwise, weights=wtss)
obama.mixed <- glm(mixed~black+hispanic+obama+(black*obama)+(hispanic*obama)+educ+lowclass+workclass+age+female+hunt_region+hunt_fund+hunt_polviews+year, family = binomial(link="logit"), data=obama.listwise, weights=wtss)
obama.struct <- glm(struct~black+hispanic+obama+(black*obama)+(hispanic*obama)+educ+lowclass+workclass+age+female+hunt_region+hunt_fund+hunt_polviews+year, family = binomial(link="logit"), data=obama.listwise, weights=wtss)


## Compare fit of Hunt specification v. Obama specification
	## Hunt's specification is like a restricted version of our Obama specification, with 3 restricted parameters
r_pcent <- 2*(as.numeric(logLik(obama.person_cent))  -  as.numeric(logLik(obama.person_cent.same)))
r_mixed <- 2*(as.numeric(logLik(obama.mixed))  -  as.numeric(logLik(obama.mixed.same)))
r_struct <- 2*(as.numeric(logLik(obama.struct))  -  as.numeric(logLik(obama.struct.same)))

Hv0_pcent <- 1-pchisq(r_pcent,df=3)
Hv0_mixed <- 1-pchisq(r_mixed,df=3)
Hv0_struct <- 1-pchisq(r_struct,df=3)

Hv0_pcent
Hv0_mixed 
Hv0_struct


####Summaries#########
summary.obama.person_cent <- summary(obama.person_cent, subset = 1:5)$coefficients
summary.obama.mixed <- summary(obama.mixed, subset = 1:5)$coefficients
summary.obama.struct <- summary(obama.struct, subset = 1:5)$coefficients


########################################################################
#########Quantities of Interest for Obama Listwise Regression###########
########################################################################

## Some of these QoI will present over time, so define:
year.obama <- c(1985,1986, 1988, 1989,1990, 1991, 1993, 1994, 1996, 1998, 2000, 2002, 2004, 2006, 2008, 2010)

## Set predictive covariates to median - model used doesn't matter since they have the same covariates
	#Allow to vary by year
Xwhite.obama <-  setx(obama.person_cent, fn=list(numeric=median), black  =  0, hispanic=0, year=year.obama)
Xblack.obama <-  setx(obama.person_cent, fn=list(numeric=median),  black  =  1, hispanic=0, year=year.obama)
Xhisp.obama  <-  setx(obama.person_cent, fn=list(numeric=median),  black  =  0, hispanic=1, year=year.obama)

## Since we cannot use sim when model not run with Zelig, compute first differences manually
	#Turn our X's into matrices
Xwhite.obama <- as.matrix(Xwhite.obama)
Xblack.obama <- as.matrix(Xblack.obama)
Xhisp.obama <- as.matrix(Xhisp.obama)

## Change value of obama to 1 for years after 2008
Xwhite.obama[,"obama"] <- as.numeric(Xblack.obama[,"year"]>=2008)
Xblack.obama[,"obama"] <- as.numeric(Xblack.obama[,"year"]>=2008)
Xhisp.obama[,"obama"] <- as.numeric(Xhisp.obama[,"year"]>=2008)

## Change value of interactions black*obama and hisp*obama
Xblack.obama[,"black:obama"] <- as.numeric(Xblack.obama[,"obama"]==1)
View(Xblack.obama)

## Predictive vectors for Blacks
Xhisp.obama[,"hispanic:obama"] <- as.numeric(Xhisp.obama[,"obama"]==1)
View(Xhisp.obama)


#################
#Person-Centered#
#################

## Simulate betas from logit results
set.seed(12345)
beta.obama.pcent <- mvrnorm(10000, mu = obama.person_cent$coef, Sigma =vcov(obama.person_cent))

##Calculate first differences for blacks v. whites
fd.obama.pcentblack <- (1 / (1 + exp(-Xblack.obama%*%t(beta.obama.pcent)))) - (1 / (1 + exp(-Xwhite.obama%*%t(beta.obama.pcent))))
meanfd.obama.pcentblack <- apply(fd.obama.pcentblack, 1, mean)
cilowfd.obama.pcentblack<- apply(fd.obama.pcentblack, 1, quantile, .025)
cihighfd.obama.pcentblack <- apply(fd.obama.pcentblack, 1, quantile, .975)


## Calculate first differences for hisp v. whites
fd.obama.pcenthisp<- (1 / (1 + exp(-Xhisp.obama%*%t(beta.obama.pcent)))) - (1 / (1 + exp(-Xwhite.obama%*%t(beta.obama.pcent))))
meanfd.obama.pcenthisp <- apply(fd.obama.pcenthisp, 1, mean)
cilowfd.obama.pcenthisp <- apply(fd.obama.pcenthisp, 1, quantile, .025)
cihighfd.obama.pcenthisp <- apply(fd.obama.pcenthisp, 1, quantile, .975)


## Plot with vertical CI for each year - Person-centered
plot(x=year.obama ,  xlab="Year", y=meanfd.obama.pcentblack, , col  =  "blue", main="Figure 2: Effect of race on probability of person-centered explanation", ylab="First differences", ylim  =  c(-0.2,0.01))
segments(x0= year.obama, x1= year.obama, y0  =  cilowfd.obama.pcentblack, y1  =  cihighfd.obama.pcentblack, col  =  "blue")
points(x=year.obama ,  y=meanfd.obama.pcenthisp,  col  =  "red")
segments(x0= year.obama, x1= year.obama, y0  =  cilowfd.obama.pcenthisp, y1  =  cihighfd.obama.pcenthisp, col  =  "red")

legend("topright", c("Blacks v. Whites","Hispanics v. Whites"),
  col=c("blue", "red"), lty=c(1,1))


#######
#Mixed#
#######

## Simulate betas from logit results
set.seed(12345)
beta.obama.mixed <- mvrnorm(10000, mu = obama.mixed$coef, Sigma =vcov(obama.mixed))

## Calculate first differences for blacks v. whites
fd.obama.mixedblack <- (1 / (1 + exp(-Xblack.obama%*%t(beta.obama.mixed)))) - (1 / (1 + exp(-Xwhite.obama%*%t(beta.obama.mixed))))
meanfd.obama.mixedblack <- apply(fd.obama.mixedblack, 1, mean)
cilowfd.obama.mixedblack<- apply(fd.obama.mixedblack, 1, quantile, .025)
cihighfd.obama.mixedblack <- apply(fd.obama.mixedblack, 1, quantile, .975)

## Calculate first differences for hisp v. whites
fd.obama.mixedhisp<- (1 / (1 + exp(-Xhisp.obama%*%t(beta.obama.mixed)))) - (1 / (1 + exp(-Xwhite.obama%*%t(beta.obama.mixed))))
meanfd.obama.mixedhisp <- apply(fd.obama.mixedhisp, 1, mean)
cilowfd.obama.mixedhisp <- apply(fd.obama.mixedhisp, 1, quantile, .025)
cihighfd.obama.mixedhisp <- apply(fd.obama.mixedhisp, 1, quantile, .975)


## Plot with vertical CI for each year - Mixed
plot(x=year.obama ,  xlab="Year", y=meanfd.obama.mixedblack, , col  =  "blue", main="Figure 2: Effect of race on probability of mixed explanation", ylab="First differences", ylim  =  c(-0.05,0.35))
segments(x0= year.obama, x1= year.obama, y0  =  cilowfd.obama.mixedblack, y1  =  cihighfd.obama.mixedblack, col  =  "blue")
points(x=year.obama ,  y=meanfd.obama.mixedhisp,  col  =  "red")
segments(x0= year.obama, x1= year.obama, y0  =  cilowfd.obama.mixedhisp, y1  =  cihighfd.obama.mixedhisp, col  =  "red")

legend("topright", c("Blacks v. Whites","Hispanics v. Whites"),
  col=c("blue", "red"), lty=c(1,1))

############
#Structural#
############

## Simulate betas from logit results
set.seed(12345)
beta.obama.struct <- mvrnorm(10000, mu = obama.struct$coef, Sigma =vcov(obama.struct))

## Calculate first differences for blacks v. whites
fd.obama.structblack <- (1 / (1 + exp(-Xblack.obama%*%t(beta.obama.struct)))) - (1 / (1 + exp(-Xwhite.obama%*%t(beta.obama.struct))))
meanfd.obama.structblack <- apply(fd.obama.structblack, 1, mean)
cilowfd.obama.structblack<- apply(fd.obama.structblack, 1, quantile, .025)
cihighfd.obama.structblack <- apply(fd.obama.structblack, 1, quantile, .975)

## Calculate first differences for hisp v. whites
fd.obama.structhisp<- (1 / (1 + exp(-Xhisp.obama%*%t(beta.obama.struct)))) - (1 / (1 + exp(-Xwhite.obama%*%t(beta.obama.struct))))
meanfd.obama.structhisp <- apply(fd.obama.structhisp, 1, mean)
cilowfd.obama.structhisp <- apply(fd.obama.structhisp, 1, quantile, .025)
cihighfd.obama.structhisp <- apply(fd.obama.structhisp, 1, quantile, .975)


## Plot with vertical CI for each year - Structural
plot(x=year.obama ,  xlab="Year", y=meanfd.obama.structblack, , col  =  "blue", main="Figure 2: Effect of race on probability of structural explanation", ylab="First differences", ylim  =  c(-0.1,0.35))
segments(x0= year.obama, x1= year.obama, y0  =  cilowfd.obama.structblack, y1  =  cihighfd.obama.structblack, col  =  "blue")
points(x=year.obama ,  y=meanfd.obama.structhisp,  col  =  "red")
segments(x0= year.obama, x1= year.obama, y0  =  cilowfd.obama.structhisp, y1  =  cihighfd.obama.structhisp, col  =  "red")

legend("topright", c("Blacks v. Whites","Hispanics v. Whites"),
  col=c("blue", "red"), lty=c(1,1))


#######################################################
########Multiple imputation for obama effect
#######################################################
set.seed(12345)
require(Amelia)
library(Zelig)

## Creating dataset for amelia that only has variables we include in the regression
amelia.obama.data <- with(obama.mat, as.data.frame(cbind(racdif_abil, racdif_motiv, racdif_educ, 
	racdif_disc, black, hispanic,educ, lowclass, workclass, age, female, hunt_region, 
	hunt_fund, hunt_polviews, year, obama)))
amelia.obama.data <- amelia.data[amelia.data$year!=1977,]


## Setting which variables are categorical
noms <- c("racdif_disc", "racdif_abil", "racdif_educ", "racdif_motiv", "hunt_region", "hunt_fund", "hunt_polviews")


## We use the amelia() function to create the five bootstrapped em imputation models. Notice the use of the cs argument and ts (not used) argument to reflect the panel nature of our data. 
## amelia.out will be a list of 5 imputed datasets, each of which can be accessed using a.out$imputations[[i]]
amelia.obama.out <- amelia(x = amelia.data, m = 5, noms = noms) # noms tells r which variables are nominal
summary(amelia.out$imputations)


## Save imputed datasets (default .csv) 
write.amelia(obj=amelia.obama.out, file.stem = "amelia.obama.out")

## Only run this section if pulling in already imputed datasets. Read in files so we don't have to run imputations every time
amelia.out <- c()
amelia.obama.out$imputations <- c()
amelia.obama.out$imputations$imp1 <- read.csv("amelia.out.obama1.csv")
amelia.obama.out$imputations$imp2 <- read.csv("amelia.out.obama2.csv")
amelia.obama.out$imputations$imp3 <- read.csv("amelia.out.obama3.csv")
amelia.obama.out$imputations$imp4 <- read.csv("amelia.out.obama4.csv")
amelia.obama.out$imputations$imp5 <- read.csv("amelia.out.obama5.csv")

## Getting rid of first column, X, which is just the row number
amelia.obama.out$imputations$imp1 <- subset(amelia.obama.out$imputations$imp1, select = c(racdif_abil, racdif_motiv, racdif_educ, 
	racdif_disc, black, hispanic, educ, lowclass, workclass, age, female, hunt_region, hunt_fund, hunt_polviews, year, obama))
amelia.obama.out$imputations$imp2 <- subset(amelia.obama.out$imputations$imp2, select = c(racdif_abil, racdif_motiv, racdif_educ, 
	racdif_disc, black, hispanic, educ, lowclass, workclass, age, female, hunt_region, hunt_fund, hunt_polviews, year, obama))
amelia.obama.out$imputations$imp3 <- subset(amelia.obama.out$imputations$imp3, select = c(racdif_abil, racdif_motiv, racdif_educ, 
	racdif_disc, black, hispanic, educ, lowclass, workclass, age, female, hunt_region, hunt_fund, hunt_polviews, year, obama))
amelia.obama.out$imputations$imp4 <- subset(amelia.obama.out$imputations$imp4, select = c(racdif_abil, racdif_motiv, racdif_educ, 
	racdif_disc, black, hispanic, educ, lowclass, workclass, age, female, hunt_region, hunt_fund, hunt_polviews, year, obama))
amelia.obama.out$imputations$imp5 <- subset(amelia.obama.out$imputations$imp5, select = c(racdif_abil, racdif_motiv, racdif_educ, 
	racdif_disc, black, hispanic, educ, lowclass, workclass, age, female, hunt_region, hunt_fund, hunt_polviews, year, obama))
## End of section where pulling in previously imputed datasets


## Defining ability
amelia.obama.out$imputations$imp1$ability <- with(amelia.obama.out$imputations$imp1, as.numeric(racdif_abil==1 & racdif_educ == 0 & racdif_disc == 0))
amelia.obama.out$imputations$imp2$ability <- with(amelia.obama.out$imputations$imp2, as.numeric(racdif_abil==1 & racdif_educ == 0 & racdif_disc == 0))
amelia.obama.out$imputations$imp3$ability <- with(amelia.obama.out$imputations$imp3, as.numeric(racdif_abil==1 & racdif_educ == 0 & racdif_disc == 0))
amelia.obama.out$imputations$imp4$ability <- with(amelia.obama.out$imputations$imp4, as.numeric(racdif_abil==1 & racdif_educ == 0 & racdif_disc == 0))
amelia.obama.out$imputations$imp5$ability <- with(amelia.obama.out$imputations$imp5, as.numeric(racdif_abil==1 & racdif_educ == 0 & racdif_disc == 0))

## Defining motivation
amelia.obama.out$imputations$imp1$motivation <- with(amelia.obama.out$imputations$imp1, as.numeric(racdif_motiv== 1 & racdif_educ == 0 & racdif_disc == 0 & racdif_abil==0))
amelia.obama.out$imputations$imp2$motivation <- with(amelia.obama.out$imputations$imp2, as.numeric(racdif_motiv== 1 & racdif_educ == 0 & racdif_disc == 0 & racdif_abil==0))
amelia.obama.out$imputations$imp3$motivation <- with(amelia.obama.out$imputations$imp3, as.numeric(racdif_motiv== 1 & racdif_educ == 0 & racdif_disc == 0 & racdif_abil==0))
amelia.obama.out$imputations$imp4$motivation <- with(amelia.obama.out$imputations$imp4, as.numeric(racdif_motiv== 1 & racdif_educ == 0 & racdif_disc == 0 & racdif_abil==0))
amelia.obama.out$imputations$imp5$motivation <- with(amelia.obama.out$imputations$imp5, as.numeric(racdif_motiv== 1 & racdif_educ == 0 & racdif_disc == 0 & racdif_abil==0))

## Defining abil_struct
amelia.obama.out$imputations$imp1$abil_struct <- with(amelia.obama.out$imputations$imp1, as.numeric(racdif_abil== 1 & (racdif_educ == 1 | racdif_disc == 1)))
amelia.obama.out$imputations$imp2$abil_struct <- with(amelia.obama.out$imputations$imp2, as.numeric(racdif_abil== 1 & (racdif_educ == 1 | racdif_disc == 1)))
amelia.obama.out$imputations$imp3$abil_struct <- with(amelia.obama.out$imputations$imp3, as.numeric(racdif_abil== 1 & (racdif_educ == 1 | racdif_disc == 1)))
amelia.obama.out$imputations$imp4$abil_struct <- with(amelia.obama.out$imputations$imp4, as.numeric(racdif_abil== 1 & (racdif_educ == 1 | racdif_disc == 1)))	
amelia.obama.out$imputations$imp5$abil_struct <- with(amelia.obama.out$imputations$imp5, as.numeric(racdif_abil== 1 & (racdif_educ == 1 | racdif_disc == 1)))

## Defining motiv_struct
amelia.obama.out$imputations$imp1$motiv_struct <- with(amelia.obama.out$imputations$imp1, as.numeric(racdif_abil== 0 & racdif_motiv== 1 & (racdif_educ == 1 | racdif_disc == 1)))
amelia.obama.out$imputations$imp2$motiv_struct <- with(amelia.obama.out$imputations$imp2, as.numeric(racdif_abil== 0 & racdif_motiv== 1 & (racdif_educ == 1 | racdif_disc == 1)))
amelia.obama.out$imputations$imp3$motiv_struct <- with(amelia.obama.out$imputations$imp3, as.numeric(racdif_abil== 0 & racdif_motiv== 1 & (racdif_educ == 1 | racdif_disc == 1)))
amelia.obama.out$imputations$imp4$motiv_struct <- with(amelia.obama.out$imputations$imp4, as.numeric(racdif_abil== 0 & racdif_motiv== 1 & (racdif_educ == 1 | racdif_disc == 1)))
amelia.obama.out$imputations$imp5$motiv_struct <- with(amelia.obama.out$imputations$imp5, as.numeric(racdif_abil== 0 & racdif_motiv== 1 & (racdif_educ == 1 | racdif_disc == 1)))

## Defining education
amelia.obama.out$imputations$imp1$education <- with(amelia.obama.out$imputations$imp1, as.numeric(racdif_abil== 0 & racdif_motiv==0 & racdif_educ == 1 & racdif_disc == 0))
amelia.obama.out$imputations$imp2$education <- with(amelia.obama.out$imputations$imp2, as.numeric(racdif_abil== 0 & racdif_motiv==0 & racdif_educ == 1 & racdif_disc == 0))
amelia.obama.out$imputations$imp3$education <- with(amelia.obama.out$imputations$imp3, as.numeric(racdif_abil== 0 & racdif_motiv==0 & racdif_educ == 1 & racdif_disc == 0))
amelia.obama.out$imputations$imp4$education <- with(amelia.obama.out$imputations$imp4, as.numeric(racdif_abil== 0 & racdif_motiv==0 & racdif_educ == 1 & racdif_disc == 0))
amelia.obama.out$imputations$imp5$education <- with(amelia.obama.out$imputations$imp5, as.numeric(racdif_abil== 0 & racdif_motiv==0 & racdif_educ == 1 & racdif_disc == 0))

## Defining discrimination
amelia.obama.out$imputations$imp1$discrimination <- with(amelia.obama.out$imputations$imp1, as.numeric(racdif_abil== 0 & racdif_motiv== 0 & racdif_disc == 1))
amelia.obama.out$imputations$imp2$discrimination <- with(amelia.obama.out$imputations$imp2, as.numeric(racdif_abil== 0 & racdif_motiv== 0 & racdif_disc == 1))
amelia.obama.out$imputations$imp3$discrimination <- with(amelia.obama.out$imputations$imp3, as.numeric(racdif_abil== 0 & racdif_motiv== 0 & racdif_disc == 1))
amelia.obama.out$imputations$imp4$discrimination <- with(amelia.obama.out$imputations$imp4, as.numeric(racdif_abil== 0 & racdif_motiv== 0 & racdif_disc == 1))
amelia.obama.out$imputations$imp5$discrimination <- with(amelia.obama.out$imputations$imp5, as.numeric(racdif_abil== 0 & racdif_motiv== 0 & racdif_disc == 1))

## Defining none
amelia.obama.out$imputations$imp1$none <- with(amelia.obama.out$imputations$imp1, as.numeric(racdif_abil == 0 & racdif_motiv == 0 & racdif_educ == 0 & racdif_disc == 0))
amelia.obama.out$imputations$imp2$none <- with(amelia.obama.out$imputations$imp2, as.numeric(racdif_abil == 0 & racdif_motiv == 0 & racdif_educ == 0 & racdif_disc == 0))
amelia.obama.out$imputations$imp3$none <- with(amelia.obama.out$imputations$imp3, as.numeric(racdif_abil == 0 & racdif_motiv == 0 & racdif_educ == 0 & racdif_disc == 0))
amelia.obama.out$imputations$imp4$none <- with(amelia.obama.out$imputations$imp4, as.numeric(racdif_abil == 0 & racdif_motiv == 0 & racdif_educ == 0 & racdif_disc == 0))
amelia.obama.out$imputations$imp5$none <- with(amelia.obama.out$imputations$imp5, as.numeric(racdif_abil == 0 & racdif_motiv == 0 & racdif_educ == 0 & racdif_disc == 0))

## Defining person_cent
amelia.obama.out$imputations$imp1$person_cent <- with(amelia.obama.out$imputations$imp1, as.numeric(ability == 1 | motivation == 1))
amelia.obama.out$imputations$imp2$person_cent <- with(amelia.obama.out$imputations$imp2, as.numeric(ability == 1 | motivation == 1))
amelia.obama.out$imputations$imp3$person_cent <- with(amelia.obama.out$imputations$imp3, as.numeric(ability == 1 | motivation == 1))
amelia.obama.out$imputations$imp4$person_cent <- with(amelia.obama.out$imputations$imp4, as.numeric(ability == 1 | motivation == 1))
amelia.obama.out$imputations$imp5$person_cent <- with(amelia.obama.out$imputations$imp5, as.numeric(ability == 1 | motivation == 1))

## Defining struct
amelia.obama.out$imputations$imp1$mixed <- with(amelia.obama.out$imputations$imp1, as.numeric(abil_struct == 1 | motiv_struct == 1))
amelia.obama.out$imputations$imp2$mixed <- with(amelia.obama.out$imputations$imp2, as.numeric(abil_struct == 1 | motiv_struct == 1))
amelia.obama.out$imputations$imp3$mixed <- with(amelia.obama.out$imputations$imp3, as.numeric(abil_struct == 1 | motiv_struct == 1))
amelia.obama.out$imputations$imp4$mixed <- with(amelia.obama.out$imputations$imp4, as.numeric(abil_struct == 1 | motiv_struct == 1))
amelia.obama.out$imputations$imp5$mixed <- with(amelia.obama.out$imputations$imp5, as.numeric(abil_struct == 1 | motiv_struct == 1))

## Defining struct
amelia.obama.out$imputations$imp1$struct <- with(amelia.obama.out$imputations$imp1, as.numeric(education == 1 | discrimination == 1))
amelia.obama.out$imputations$imp2$struct <- with(amelia.obama.out$imputations$imp2, as.numeric(education == 1 | discrimination == 1))
amelia.obama.out$imputations$imp3$struct <- with(amelia.obama.out$imputations$imp3, as.numeric(education == 1 | discrimination == 1))
amelia.obama.out$imputations$imp4$struct <- with(amelia.obama.out$imputations$imp4, as.numeric(education == 1 | discrimination == 1))
amelia.obama.out$imputations$imp5$struct <- with(amelia.obama.out$imputations$imp5, as.numeric(education == 1 | discrimination == 1))


## Run model again using zelig since it can handle multiply imputed datasets
imp.obama.ability <- zelig(ability~black+hispanic+obama+(black*obama)+(hispanic*obama)+educ+lowclass+workclass+age+female+hunt_region+hunt_fund+hunt_polviews+year, family = binomial(link="logit"), data=amelia.obama.out$imputations, model = "logit")
imp.obama.motivation <- zelig(motivation~black+hispanic+obama+(black*obama)+(hispanic*obama)+educ+lowclass+workclass+age+female+hunt_region+hunt_fund+hunt_polviews+year, family = binomial(link="logit"), data=amelia.obama.out$imputations, model = "logit")
imp.obama.abil_struct <- zelig(abil_struct~black+hispanic+obama+(black*obama)+(hispanic*obama)+educ+lowclass+workclass+age+female+hunt_region+hunt_fund+hunt_polviews+year, family = binomial(link="logit"), data=amelia.obama.out$imputations, model = "logit")
imp.obama.motiv_struct <- zelig(motiv_struct~black+hispanic+obama+(black*obama)+(hispanic*obama)+educ+lowclass+workclass+age+female+hunt_region+hunt_fund+hunt_polviews+year, family = binomial(link="logit"), data=amelia.obama.out$imputations, model = "logit")
imp.obama.education <- zelig(education~black+hispanic+obama+(black*obama)+(hispanic*obama)+educ+lowclass+workclass+age+female+hunt_region+hunt_fund+hunt_polviews+year, family = binomial(link="logit"), data=amelia.obama.out$imputations, model = "logit")
imp.obama.discrimination <- zelig(discrimination~black+hispanic+obama+(black*obama)+(hispanic*obama)+educ+lowclass+workclass+age+female+hunt_region+hunt_fund+hunt_polviews+year, family = binomial(link="logit"), data=amelia.obama.out$imputations, model = "logit")
imp.obama.none <- zelig(none~black+hispanic+obama+(black*obama)+(hispanic*obama)+educ+lowclass+workclass+age+female+hunt_region+hunt_fund+hunt_polviews+year, family = binomial(link="logit"), data=amelia.obama.out$imputations, model = "logit")
imp.obama.person_cent <- zelig(person_cent~black+hispanic+obama+(black*obama)+(hispanic*obama)+educ+lowclass+workclass+age+female+hunt_region+hunt_fund+hunt_polviews+year, family = binomial(link="logit"), data=amelia.obama.out$imputations, model = "logit")
imp.obama.mixed <- zelig(mixed~black+hispanic+obama+(black*obama)+(hispanic*obama)+educ+lowclass+workclass+age+female+hunt_region+hunt_fund+hunt_polviews+year, family = binomial(link="logit"), data=amelia.obama.out$imputations, model = "logit")
imp.obama.struct <- zelig(struct~black+hispanic+obama+(black*obama)+(hispanic*obama)+educ+lowclass+workclass+age+female+hunt_region+hunt_fund+hunt_polviews+year, family = binomial(link="logit"), data=amelia.obama.out$imputations, model = "logit")


## The missingness map gives an overall sense of the shape and extent of the missingness.
missmap(amelia.obama.out)

## Generating table to compare to hunt's results
imp.obama.ability <- summary(imp.obama.ability, subset = 1:5)$coefficients
imp.obama.motivation <- summary(imp.obama.motivation, subset = 1:5)$coefficients
imp.obama.abil_struct <- summary(imp.obama.abil_struct, subset = 1:5)$coefficients
imp.obama.motiv_struct <- summary(imp.obama.motiv_struct, subset = 1:5)$coefficients
imp.obama.education <- summary(imp.obama.education, subset = 1:5)$coefficients
imp.obama.discrimination <- summary(imp.obama.discrimination, subset = 1:5)$coefficients
imp.obama.none <- summary(imp.obama.none, subset = 1:5)$coefficients
imp.obama.person_cent <- summary(imp.obama.person_cent, subset = 1:5)$coefficients
imp.obama.mixed <- summary(imp.obama.mixed, subset = 1:5)$coefficients
imp.obama.struct <- summary(imp.obama.struct, subset = 1:5)$coefficients

table.ability <- write.csv(cbind(imp.obama.ability,"", summary.obama.ability), 
	"Obama imputed versus listwise - Ability.csv",row.names=T)
table.motivation <- write.csv(cbind(imp.obama.motivation,"", summary.obama.motivation), 
	"Obama imputed versus listwise - motiv.csv",row.names=T)
table.abil_struct<- write.csv(cbind(imp.obama.abil_struct,"", summary.obama.abil_struct), 
	"Obama imputed versus listwise - abil_struct.csv",row.names=T)
table.motiv_struct<- write.csv(cbind(imp.obama.motiv_struct,"", summary.obama.motiv_struct),
	"Obama imputed versus listwise - motiv_struct.csv",row.names=T)
table.education <- write.csv(cbind(imp.obama.education ,"", summary.obama.education), 
	"Obama imputed versus listwise - education.csv",row.names=T)
table.discrimination <- write.csv(cbind(imp.obama.discrimination, "", summary.obama.discrimination), 
	"Obama imputed versus listwise - discrimination.csv",row.names=T)
table.none <- write.csv(cbind(imp.obama.none, "", summary.obama.none), 
	"Obama imputed versus listwise - none.csv",row.names=T)
table.person_cent <- write.csv(cbind(imp.obama.person_cent, "", summary.obama.person_cent), 
	"Obama imputed versus listwise - person_cent.csv",row.names=T)
table.mixed <- write.csv(cbind(imp.obama.mixed, "", summary.obama.mixed), 
	"Obama imputed versus listwise - mixed.csv",row.names=T)
table.struct <- write.csv(cbind(imp.obama.struct, "", summary.obama.struct), 
	"Obama imputed versus listwise - struct.csv",row.names=T)

